有限差分法求解偏微分方程MATLAB

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

2、有限差分法的基本思想是把連續(xù)的定解區(qū)域用有限個(gè)離散點(diǎn)構(gòu)成的網(wǎng)格來(lái)代替；把連續(xù)定解區(qū)域上的連續(xù)變量的函數(shù)用在網(wǎng)格上定義的離散變量函數(shù)來(lái)近似；把原方程和定解條件中的微商用差商來(lái)近似，積分用積分和來(lái)近似，于是原微分方程和定解條件就近似地代之以代數(shù)方程組，即有限差分方程組，解此方程組就可以得到原問(wèn)題在離散點(diǎn)上的近似解。推導(dǎo)差分方程的過(guò)程中需要用到的泰勒展開(kāi)公式如下： (2-1)求解區(qū)域的網(wǎng)格劃分步長(zhǎng)參數(shù)如下： (2-2)2.1 古典顯格式2.1.1 古典顯格式的推導(dǎo)由泰勒展開(kāi)公式將對(duì)時(shí)間展開(kāi)得 (2-3)當(dāng)時(shí)有 (2-4)得到對(duì)時(shí)間的一階偏導(dǎo)數(shù) (2-5)由泰勒展開(kāi)公式將對(duì)位置展開(kāi)得 (2-6)當(dāng)時(shí)，

3、代入式(2-6)得 (2-7)因?yàn)椋肷鲜降?(2-8)得到對(duì)位置的二階偏導(dǎo)數(shù) (2-9)將式(2-5)、(2-9)代入一般形式的拋物線型偏微分方程得(2-10)為了方便我們可以將式(2-10)寫(xiě)成 (2-11) (2-12)最后得到古典顯格式的差分格式為 (2-13)，古典顯格式的差分格式的截?cái)嗾`差是。 古典顯格式穩(wěn)定性分析古典顯格式(2-13)寫(xiě)成矩陣形式為 (2-14)上面的C矩陣的特征值是： (2-15)使，即結(jié)論：當(dāng)時(shí)，所以古典顯格式是穩(wěn)定的。2.2 古典隱格式2.2.1 古典隱格式的推導(dǎo)將代入式 (2-3)得 (2-16) (2-17)得到對(duì)時(shí)間的一階偏導(dǎo)數(shù) (2-18)將式(2

4、-9)、(2-18)原方程得到(2-19)為了方便把(2-19)寫(xiě)成 (2-20) (2-21)最后得到古典隱格式的差分格式為 (2-22)，古典隱格式的差分格式的截?cái)嗾`差是。 古典隱格式穩(wěn)定性分析將古典隱格式(2-22)寫(xiě)成矩陣形式如下 (2-23)誤差傳播方程 (2-24)所以誤差方程的系數(shù)矩陣為使，顯然 恒成立。結(jié)論：對(duì)于，即任意網(wǎng)格比下，古典隱格式是絕對(duì)穩(wěn)定的。2.3 Richardson格式 Richardson格式的推導(dǎo)將，代入式(2-3)得 (2-25)即 (2-26)由此得到可得 (2-27)將式(2-9) 、(2-27)代入原方程得到下式 (2-28)為了方便可以把式(2-2

5、8)寫(xiě)成 (2-29)即 (2-30)最后得到Richardson顯格式的差分格式為 (2-31)，古典顯格式的差分格式的截?cái)嗾`差是。 Richardson穩(wěn)定性分析將Richardson顯格式(2-31)寫(xiě)成如下矩陣形式 (2-32)誤差傳播方程矩陣形式 (2-33)再將上面的方程組寫(xiě)成矩陣形式 (2-34)系數(shù)矩陣的特征值是 (2-35)解得特征值為 (2-36) (恒成立) (2-37)結(jié)論：上式對(duì)任意的網(wǎng)比都恒成立，即Richardson格式是絕對(duì)不穩(wěn)定的。4. Crank-Nicholson格式 Crank-Nicholson格式的推導(dǎo)將代入式(2-9)得 (2-40)即 (2-41

6、)得到如下方程 (2-42)所以處的一階偏導(dǎo)數(shù)可以用下式表示： (2-43)將,代入式(2-6)可以得到式(2-9)；同理,代入式(2-6)可以得到 (2-44)所以，處的二階偏導(dǎo)數(shù)用式(2-6)、(2-44)表示： (2-45)所以，處的函數(shù)值可用下式表示： (2-46)原方程變?yōu)椋?(2-47)將差分格式代入上式得： (2-48)為了方便寫(xiě)成： (2-49)最后得到Crank-Nicholson格式的差分格式為 (2-50)，Crank-Nicholson格式的差分格式的截?cái)嗾`差是。 Crank-Nicholson穩(wěn)定性分析Crank-Nicholson格式寫(xiě)成矩陣形式如下： (2-51)

7、誤差傳播方程是： (2-52)可以得到： (2-53) (2-54)使 即 (2-55) (2-56) (2-57) (2-58)上式恒成立。結(jié)論：Crank-Nicholson格式對(duì)任意網(wǎng)格比也是絕對(duì)穩(wěn)定的。5. Du Fort Frankle格式（Richardson格式的改進(jìn)）將代入式(2-31)并化簡(jiǎn)得到Du Fort Frankle： (2-59) (2-60)可以證明Du Fort Frankle格式是絕對(duì)穩(wěn)定的。因?yàn)榇烁袷绞荝ichardson格式的改進(jìn)格式，因此截?cái)嗾`差還是。3.6 總結(jié)(1) 古典顯格式古典顯格式的差分格式為截?cái)嗾`差：。穩(wěn)定性：當(dāng)時(shí)，古典顯格式穩(wěn)定。(2) 古

8、典隱格式古典隱格式的差分格式為截?cái)嗾`差：。穩(wěn)定性：任意網(wǎng)格比古典隱格式絕對(duì)穩(wěn)定。(3) Richardson顯格式Richardson顯格式的差分格式為截?cái)嗾`差：。穩(wěn)定性：任意網(wǎng)格比Richardson格式絕對(duì)不穩(wěn)定。(4) Crank-Nicholson格式Crank-Nicholson格式的差分格式為截?cái)嗾`差：。穩(wěn)定性：Crank-Nicholson格式對(duì)任意網(wǎng)格比絕對(duì)穩(wěn)定。(5) Du Fort Frankle格式截?cái)嗾`差：。穩(wěn)定性：Du Fort Frankle格式對(duì)任意網(wǎng)格比絕對(duì)穩(wěn)定。三、MATLAB實(shí)現(xiàn)五種差分格式對(duì)偏微分方程的求解及誤差分析3.1 精確數(shù)值解上述偏微分方程的精確解

9、是區(qū)域取值范圍：。用MATLAB對(duì)精確解進(jìn)行編程畫(huà)出三維圖像精確解程序如下：close allclcx,t=meshgrid(0:0.01:1,0:0.001:0.2)u=exp(-pi*pi*t).*sin(pi*x)mesh(x,t,u);surf(x,t,u);xlabel(x),ylabel(t),zlabel(u);title(精確數(shù)值解);shading interp圖3-1 精確數(shù)值解的三維圖(a) 精確數(shù)值解X-Y平面 (b) 精確數(shù)值解X-Z平面(c) 精確數(shù)值解Y-Z平面圖3-2 精確數(shù)值解的三個(gè)平面圖3.2 五種差分格式MATLAB程序3.2.1 古典顯格式close a

10、llclcT=0.2X=1.0M=41N=11u=zeros(M,N); %構(gòu)造一個(gè)M行N列的矩陣用于存放時(shí)間t和變量xra=(T/(M-1)/(X/(N-1)2); %網(wǎng)格比f(wàn)printf(穩(wěn)定性系數(shù) S=ra 為:n);disp(ra); % 顯示網(wǎng)格比f(wàn)or i=2:N-1 xx=(i-1)*(X/(N-1); u(1,i)=sin(pi*xx); end; %即t=0時(shí)刻賦值，邊界條件for k=1:M u(k,1)=0; u(k,N)=0;end; % x=0,x=1處的邊界條件for k=1:M-1 %矩陣是從y軸表示行k，x軸表示列i。由行開(kāi)始 for i=2:N-1 u(k+1

11、,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); %將區(qū)域劃分成網(wǎng)格對(duì)每個(gè)點(diǎn)賦值再畫(huà)圖surf(x,t,u); xlabel(x),ylabel(t),zlabel(u);title( 古典顯格式); %此程序得到的是圖3-3圖3-3古典顯格式程序結(jié)果圖（r=0.5）圖3-4精確數(shù)值解、古典顯格式程序結(jié)果的Y-Z平面圖（r=0.5）圖3-5古典顯格式在取不同網(wǎng)格比時(shí)的誤差傳播結(jié)果圖圖3-6 不同時(shí)間取值時(shí)精確解、與古典顯格式的值對(duì)比圖（網(wǎng)

12、格比r=0.5）紅線表示精確解、藍(lán)色線表示差分格式的解圖3-1、圖3-3對(duì)比可以看出，精確解和古典顯格式（網(wǎng)格比r=0.5）的圖形是一致的。圖3-4精確數(shù)值解、古典顯格式的Y-Z平面圖結(jié)果可以看出古典顯格式的邊界值和精確解一樣。圖3-5是r分別取0.245、0.5、0.72、1.125時(shí)的誤差傳播圖像，邊界位置網(wǎng)格數(shù)為5處的誤差為0.01得到的，可以看出r小于等于0.5是穩(wěn)定的；但是r大于0.5出現(xiàn)不穩(wěn)定現(xiàn)象。很好的驗(yàn)證了古典顯格式穩(wěn)定性。3.2.2 古典隱格式close allclcT=0.2X=1.0M=41N=21ra=(T/(M-1)/(X/(N-1)2); %網(wǎng)格比f(wàn)printf(穩(wěn)

13、定性系數(shù) S=ra 為:n);disp(ra); %顯示網(wǎng)格比u=zeros(M,N); %構(gòu)造一個(gè)M行N列的矩陣用于存放時(shí)間t和變量x for i=2:N-1 xx=(i-1)*(X/(N-1); u(1,i)=sin(pi*xx); %t=0時(shí)刻的賦值，邊界條件 end; for k=1:M u(k,1)=0; u(k,N)=0; end; % x=0,x=1處的邊界條件A=zeros(N-1,N-1); %隱格式的矩陣形式中的A矩陣賦值 A(1,1)=1+2*ra; A(N-1,N-1)=1+2*ra; A(1,2)=-ra; A(N-1,N-2)=-ra; for m=2:N-2 A(

14、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);for i=1:N-1 d(i,1)=sin(pi*(i-1)*(1.0/(N-1); end %隱格式右邊矩陣賦值%以下循環(huán)對(duì)矩陣A進(jìn)行LU分解U(1,1)=A(1,1); for i=2:n L(i,i-1)=A(i,i-1)/U(i-1,i-1); U(i-1,i)=A(i-1,i);%U的上次對(duì)角線即為A的上次對(duì)

15、角線 U(i,i)=A(i,i)-L(i,i-1)*U(i-1,i); endfor k=1:M-1 %外層大循環(huán)是求時(shí)間網(wǎng)格2,3，M的求解u%以下是追趕法的求解過(guò)程%-追的過(guò)程-即Ly=d的求解yy(1,1)=d(1,1);for i=2:n y(i,1)=d(i,1)-L(i,i-1)*y(i-1,1);end%-趕的過(guò)程-即Ux=y的求解xx(n,1)=y(n,1)/U(n,n);for i=n-1:-1:1 x(i,1)=(y(i,1)-U(i,i+1)*x(i+1,1)/U(i,i);end %追趕法結(jié)束 for i=1:n u(k+1,i)=x(i,1) end d=zeros(

16、N-1,1); %更新右邊矩陣 d=x %每次外循環(huán)更換右邊矩陣end for k=1:M u(k,1)=0; end;disp(u); % 顯示差分法求得的值 x,t=meshgrid(1:N,1:M); %將區(qū)域劃分成網(wǎng)格對(duì)每個(gè)點(diǎn)賦值再畫(huà)圖surf(x,t,u);xlabel(x),ylabel(t),zlabel(u);title(古典隱格式); %此程序得到圖是3-7圖3-7古典隱格式程序結(jié)果圖（r=2）圖3-8精確數(shù)值解、古典隱格式程序結(jié)果的Y-Z平面圖（r=2）圖3-9古典隱格式在取不同網(wǎng)格比時(shí)的誤差傳播結(jié)果圖圖3-10 不同時(shí)間取值時(shí)精確解、與古典隱格式的值對(duì)比圖（網(wǎng)格比r=2）

17、紅線表示精確解、藍(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時(shí)的誤差傳播圖像，邊界位置網(wǎng)格數(shù)為5處的誤差為0.01得到的，可以看出r取不同的值時(shí)都是穩(wěn)定的；即古典隱格式對(duì)任意的網(wǎng)格比穩(wěn)定性。從圖3-10可以看出隱格式隨著時(shí)間的增加，差分格式計(jì)算的結(jié)果和精確結(jié)果越來(lái)越大；隱格式雖然對(duì)任意網(wǎng)格比都是穩(wěn)定的，但是計(jì)算的精確度是它的不足。3.2.3 Richardson顯格式close allclcT=0

18、.2X=1.0000M=41N=11ra=(T/(M-1)/(X/(N-1)2);fprintf(穩(wěn)定性系數(shù) S=ra 為:n);disp(ra);u=zeros(M,N); %構(gòu)造一個(gè)M行N列的矩陣 for i=2:N-1 xx=(i-1)*(X/(N-1); u(1,i)=sin(pi*xx); %邊界條件 end; for k=1:M u(k,1)=0; u(k,N)=0; % 邊界條件 end;%因?yàn)镽ichadson格式需要知道前兩行的值，第二行值我們采用隱格式求得%下面是隱格式求解第二行，和隱格式的程序一樣，只是求一行，此處不再做注釋A=zeros(N-1,N-1); A(1,1)

19、=1+2*ra; A(N-1,N-1)=1+2*ra; A(1,2)=-ra; A(N-1,N-2)=-ra; for m=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); for i=2:n L(i,i-1)=A(i,i-1)/U(i-1,i-1); U(i-1,i)=A(i-1,i); %U的上次對(duì)角線即為A的上次對(duì)角線 U(i,i)=A(i,i)-L(i,i-1)*U(i-1,i); endy(1,1

20、)=u(1,1);for i=2:n y(i,1)=u(1,i)-L(i,i-1)*y(i-1,1);endx(n,1)=y(n,1)/U(n,n);for i=n-1:-1:1 x(i,1)=(y(i,1)-U(i,i+1)*x(i+1,1)/U(i,i);end u(2,1)=0;for i=2:N-1 u(2,i)=x(i,1)end%通過(guò)隱格式已求得第二行的值u(2,i)，下面是Richardson格式的過(guò)程for k=2:M-1 %時(shí)間軸第3行開(kāi)始到第M行 for i=2:N-1 %位置網(wǎng)格點(diǎn) u(k+1,i)=2*ra*(u(k,i-1)-2*u(k,i)+u(k,i+1)+u(k

21、-1,i); %Richardson格式 endenddisp(u); %顯示求解的值x,t=meshgrid(1:N,1:M); %區(qū)域劃分進(jìn)行賦值畫(huà)圖surf(x,t,u);xlabel(x),ylabel(t),zlabel(u);title(Richardson格式); %此程序得到的圖形是圖3-11圖3-11 Richardson顯格式程序結(jié)果圖（r=0.5）圖3-12精確數(shù)值解、Richardson顯格式程序結(jié)果的Y-Z平面圖（r=0.5）圖3-13 Richardson顯格式在取不同網(wǎng)格比時(shí)都不穩(wěn)定的結(jié)果圖圖3-11是 Richardson顯格式在r=0.5時(shí)的程序結(jié)果圖，可以看

22、出不穩(wěn)定。圖3-12是精確數(shù)值解、Richardson顯格式程序結(jié)果的Y-Z平面圖。從圖3-13可以看出Richardson顯格式在取不同網(wǎng)格比時(shí)都出現(xiàn)不穩(wěn)定現(xiàn)象，和理論結(jié)果相一致。所以說(shuō)Richardson顯格式絕對(duì)不穩(wěn)定，這種差分格式不能使用。后面有改進(jìn)的格式D-F格式。3.2.4 Crank-Nicholson格式close allclcT=0.2X=1.0000M=41N=21ra=(T/(M-1)/(X/(N-1)2);fprintf(穩(wěn)定性系數(shù) S=ra 為:n);disp(ra);u=zeros(M,N); %構(gòu)造一個(gè)M行N列的矩陣用于存放時(shí)間t和變量x%disp(u);for

23、i=2:N-1 xx=(i-1)*(X/(N-1); u(1,i)=sin(pi*xx); % 邊界條件 end; for k=1:M u(k,1)=0; u(k,N)=0; % 邊界條件 end; %Crank-Nicholson格式需要兩個(gè)矩陣，下面的A、B，參照Crank-Nicholson格式的矩陣形式A=zeros(N-1,N-1); %下面對(duì)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; for m=2:N-2 A(m,m-1)=-ra/2; A(m,m)=1+ra; A(m,m+1)=-r

24、a/2;end;B=zeros(N-1,N-1); %下面對(duì)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;for m=2:N-2 B(m,m-1)=ra/2; B(m,m)=1-ra; B(m,m+1)=ra/2;end;%以下是追趕法求u值d=zeros(N-1,1); %首先填充右邊向量，然后進(jìn)行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); for i=2:n L(i,i-1)=A(i,i-

25、1)/U(i-1,i-1); U(i-1,i)=A(i-1,i); %U的上次對(duì)角線即為A的上次對(duì)角線 U(i,i)=A(i,i)-L(i,i-1)*U(i-1,i); endfor i=1:N-1 d(i,1)=sin(pi*(i-1)*(1.0/(N-1); endfor k=1:M-1 %按層外圍大循環(huán)，即時(shí)間步長(zhǎng) 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); for i=2:N-2 m(i,1)=B(i,i-1)*d(i-1,1)+B

26、(i,i)*d(i,1)+B(i,i+1)*d(i+1,1); end %以上是右邊矩陣的填充更新%-追-求解Ly=b，其中b是原來(lái)的列向量矩陣乘上B系數(shù)矩陣得到的y(1,1)=m(1,1);for i=2:n 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);for i=n-1:-1:1 x(i,1)=(y(i,1)-U(i,i+1)*x(i+1,1)/U(i,i);end for i=1:n u(k+1,i)=x(i,1) end d=zeros(N-1,1); %右邊向量 d=xend for k=1:M

27、u(k,1)=0; end;disp(u); % u的值全部求出，以下畫(huà)圖x,t=meshgrid(1:N,1:M);surf(x,t,u);xlabel(x),ylabel(t),zlabel(u);title( Crank-Nicholson格式); %此程序得到的圖像是圖3-14圖3-14 Crank-Nicholson格式程序結(jié)果圖（r=2）圖3-15精確數(shù)值解、Crank-Nicholson格式程序結(jié)果的Y-Z平面圖（r=2）圖3-16 Crank-Nicholson格式在取不同網(wǎng)格比時(shí)的誤差傳播結(jié)果圖圖3-17 不同時(shí)間取值時(shí)精確解、與C-N格式的解對(duì)比圖（網(wǎng)格比r=2）紅線表示精

28、確解、藍(lán)色線表示差分格式的解圖3-14是程序運(yùn)行得到的Crank-Nicholson格式在網(wǎng)格比取r=2的結(jié)果，和精確解圖像一致。在t=0.2時(shí)從圖3-15 的Y-Z面可以看出結(jié)果還是有一定的誤差。理論上Crank-Nicholson格式對(duì)任意的網(wǎng)格比也是穩(wěn)定的，從圖3-16可以看出在r取0.245、0.5、0.72、1.125誤差傳播圖像可以看出誤差不會(huì)擴(kuò)撒。圖3-17是不同時(shí)間取值時(shí)精確解、與C-N格式r=2時(shí)的解對(duì)比圖，計(jì)算還是具有誤差。 Du Fort Frankle格式close allclcT=0.2X=1.0000M=41N=21ra=(T/(M-1)/(X/(N-1)2);fp

29、rintf(穩(wěn)定性系數(shù) S=ra 為:n);disp(ra);u=zeros(M,N); %構(gòu)造一個(gè)M行N列的矩陣 for i=2:N-1 xx=(i-1)*(X/(N-1); u(1,i)=sin(pi*xx); %i表示x軸，k表示y軸(即t) end; for k=1:M %其實(shí)初始矩陣已經(jīng)將i=1和i=N列的初值賦零了 u(k,1)=0; u(k,N)=0; end;%第二行用隱格式求得A=zeros(N-1,N-1); %下面對(duì)A進(jìn)行填充賦值 A(1,1)=1+2*ra; A(N-1,N-1)=1+2*ra; A(1,2)=-ra; %第一行第二個(gè)和最后一行倒數(shù)第二個(gè)一個(gè)賦值 A(N

30、-1,N-2)=-ra; for m=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); for i=2:n L(i,i-1)=A(i,i-1)/U(i-1,i-1); U(i-1,i)=A(i-1,i); %U的上次對(duì)角線即為A的上次對(duì)角線 U(i,i)=A(i,i)-L(i,i-1)*U(i-1,i); end%-追-y(1,1)=u(1,1);for i=2:n y(i,1)=u(1,i)-L(i,i-

31、1)*y(i-1,1);end%-趕-x(n,1)=y(n,1)/U(n,n);for i=n-1:-1:1 x(i,1)=(y(i,1)-U(i,i+1)*x(i+1,1)/U(i,i);end u(2,1)=0;for i=2:N-1 u(2,i)=x(i,1)end %第二行求出，下面用D-F差分格式for k=2:M-1 for i=2:N-1 u(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(1:N,1:M);surf(x,t,u);xlabel(x

32、),ylabel(t),zlabel(u);title(Du Fort Frankle格式); %此程序?yàn)镽ichardson格式的改進(jìn)，得到圖3-18圖3-18 Du Fort Frankle格式程序結(jié)果圖（r=2）圖3-19精確數(shù)值解、Du Fort Frankle格式程序結(jié)果的Y-Z平面圖（r=2）圖3-20 Du Fort Frankle格式在取不同網(wǎng)格比時(shí)的誤差傳播結(jié)果圖圖3-21 不同時(shí)間取值時(shí)精確解、與D-F格式的解對(duì)比圖（網(wǎng)格比r=2）紅線表示精確解、藍(lán)色線表示差分格式的解D-F格式是Richardson格式（絕對(duì)不穩(wěn)定）的改進(jìn)，從圖3-18可以看出當(dāng)r時(shí)D-F格式是穩(wěn)定的；圖

33、3-20表示D-F格式網(wǎng)格比r為0.245、0.5、0.72、1.125時(shí)誤差傳播圖像，不同的網(wǎng)格比誤差都不會(huì)擴(kuò)撒。說(shuō)明這種格式是穩(wěn)定的，和理論上的結(jié)果相一致。圖3-21是不同時(shí)間取值時(shí)精確解、與D-F格式的解對(duì)比，隨時(shí)間的增加計(jì)算的值和差分得到的值有誤差。此種格式雖然是絕對(duì)穩(wěn)定的，但是計(jì)算的精度還是有待提升。四、結(jié)論及感想從程序得出的結(jié)果與精確解的對(duì)比來(lái)看和理論是上的結(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格式理論上是絕對(duì)不穩(wěn)定的，從編程的結(jié)

34、果在取不同的網(wǎng)格比Richardon格式都是不穩(wěn)定的，理論和結(jié)果一致。理論上對(duì)不穩(wěn)定格式進(jìn)行改進(jìn)使其穩(wěn)定，比如得到的D-F格式，取不同的格式網(wǎng)格比都是穩(wěn)定的，很好的驗(yàn)證了改進(jìn)的格式的穩(wěn)定性。本次報(bào)告首先推導(dǎo)了一維拋物線型偏微分方程的一般差分格式。分別是古典顯格式、古典隱格式、Richardson顯格式、C-N格式進(jìn)版的Richardson格式D-F格式。推導(dǎo)中得到各種格式的截?cái)嗾`差，從理論上分析了各種格式的穩(wěn)定性。對(duì)于不穩(wěn)定的格式進(jìn)行修改得到穩(wěn)定的格式，即Richardson格式的修改。通過(guò)MATLAB編程實(shí)現(xiàn)了各種格式的程序?qū)崿F(xiàn)。用實(shí)驗(yàn)的方法來(lái)驗(yàn)證理論結(jié)果，能更好的理解各種差分格式的穩(wěn)定性及

35、操作過(guò)程。報(bào)告中的程序都是基本程序，誤差圖與二維x-u的圖像（見(jiàn)附錄）都是在基本程序的基礎(chǔ)上對(duì)參數(shù)的修改得到的圖像。通過(guò)這次報(bào)告的完成學(xué)到了很多的東西。首先，對(duì)編程有了進(jìn)一步的了解；尤其是使用MATLAB編程。在這個(gè)過(guò)程中也遇到了很多的問(wèn)題，通過(guò)查閱資料并利用網(wǎng)絡(luò)資源尋求解決辦法。其次，對(duì)于差分格式在程序上的實(shí)現(xiàn)并不是簡(jiǎn)單的書(shū)寫(xiě)，需要步步銜接好；比如好幾種格式都用到差分格式的矩陣形式，尤其是隱式格式不能直接求解，需要應(yīng)用追趕法進(jìn)行求解；編寫(xiě)過(guò)程中通過(guò)直接求解得出的結(jié)果都不正確，通過(guò)追趕法才能得到正確的結(jié)果。最后，我們專(zhuān)業(yè)是電磁場(chǎng)與微波技術(shù)且偏計(jì)算，經(jīng)常遇到偏微分方程的求解；所以這次試驗(yàn)毫無(wú)疑問(wèn)

36、的對(duì)我們理解有限差分法和求解方程提供了很大的幫助。最后，感謝張老師在課堂上的知識(shí)的講解；同時(shí)也感謝在完成報(bào)告期間對(duì)我提供幫助的同學(xué)！附 錄誤差傳播三維圖（以顯格式（圖3-5）為例）：close allclcT=0.2X=1.0M=41N=8u=zeros(M,N); %構(gòu)造一個(gè)M行N列的矩陣用于存放時(shí)間t和變量xra=(T/(M-1)/(X/(N-1)2); %網(wǎng)格比f(wàn)printf(穩(wěn)定性系數(shù) S=ra 為:n);disp(ra); % 顯示網(wǎng)格比u(1,5)=0.01;%即t=0時(shí)刻賦值，誤差0.01for k=1:M u(k,1)=0; u(k,N)=0;end; % x=0,x=1處的邊

37、界條件for k=1:M-1 %矩陣是從y軸表示行k，x軸表示列i。由行開(kāi)始 for i=2:N-1 u(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); %將區(qū)域劃分成網(wǎng)格對(duì)每個(gè)點(diǎn)賦值再畫(huà)圖subplot(2,2,1)surf(x,t,u); xlabel(x),ylabel(t),zlabel(u);title(古典顯格式-誤差傳播 r=0.245);clcT=0.2X=1.0M=41N=11u=zeros(M,N);ra=(T

38、/(M-1)/(X/(N-1)2);fprintf(穩(wěn)定性系數(shù) S=ra 為:n);disp(ra); u(1,5)=0.01; %即t=0時(shí)刻賦值，誤差0.01 for k=1:M u(k,1)=0; u(k,N)=0;end;for k=1:M-1 for i=2:N-1 u(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(古典顯格式-誤差傳播

39、r=0.5);clcT=0.2X=1.0M=41N=13u=zeros(M,N);ra=(T/(M-1)/(X/(N-1)2);fprintf(穩(wěn)定性系數(shù) S=ra 為:n);disp(ra);u(1,5)=0.01; %即t=0時(shí)刻賦值，誤差0.01for k=1:M u(k,1)=0; u(k,N)=0;end;for k=1:M-1 for i=2:N-1 u(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,3)surf(x,t,u); xlabel(x