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

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

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

3、代入式(2-6)得 (2-7)因為，代入上式得 (2-8)得到對位置的二階偏導數(shù) (2-9)將式(2-5)、(2-9)代入一般形式的拋物線型偏微分方程得(2-10)為了方便我們可以將式(2-10)寫成 (2-11) (2-12)最后得到古典顯格式的差分格式為 (2-13)，古典顯格式的差分格式的截斷誤差是。 古典顯格式穩(wěn)定性分析古典顯格式(2-13)寫成矩陣形式為 (2-14)上面的C矩陣的特征值是： (2-15)使，即結論：當時，所以古典顯格式是穩(wěn)定的。2.2 古典隱格式2.2.1 古典隱格式的推導將代入式 (2-3)得 (2-16) (2-17)得到對時間的一階偏導數(shù) (2-18)將式(2

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

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

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

7、誤差傳播方程是： (2-52)可以得到： (2-53) (2-54)使 即 (2-55) (2-56) (2-57) (2-58)上式恒成立。結論：Crank-Nicholson格式對任意網(wǎng)格比也是絕對穩(wěn)定的。5. Du Fort Frankle格式（Richardson格式的改進）將代入式(2-31)并化簡得到Du Fort Frankle： (2-59) (2-60)可以證明Du Fort Frankle格式是絕對穩(wěn)定的。因為此格式是Richardson格式的改進格式，因此截斷誤差還是。3.6 總結(1) 古典顯格式古典顯格式的差分格式為截斷誤差：。穩(wěn)定性：當時，古典顯格式穩(wěn)定。(2) 古

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

9、是區(qū)域取值范圍：。用MATLAB對精確解進行編程畫出三維圖像精確解程序如下：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ù)值解的三個平面圖3.2 五種差分格式MATLAB程序3.2.1 古典顯格式close a

10、llclcT=0.2X=1.0M=41N=11u=zeros(M,N); %構造一個M行N列的矩陣用于存放時間t和變量xra=(T/(M-1)/(X/(N-1)2); %網(wǎng)格比fprintf(穩(wěn)定性系數(shù) S=ra 為:n);disp(ra); % 顯示網(wǎng)格比for i=2:N-1 xx=(i-1)*(X/(N-1); u(1,i)=sin(pi*xx); end; %即t=0時刻賦值，邊界條件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。由行開始 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)格對每個點賦值再畫圖surf(x,t,u); xlabel(x),ylabel(t),zlabel(u);title( 古典顯格式); %此程序得到的是圖3-3圖3-3古典顯格式程序結果圖（r=0.5）圖3-4精確數(shù)值解、古典顯格式程序結果的Y-Z平面圖（r=0.5）圖3-5古典顯格式在取不同網(wǎng)格比時的誤差傳播結果圖圖3-6 不同時間取值時精確解、與古典顯格式的值對比圖（網(wǎng)

12、格比r=0.5）紅線表示精確解、藍色線表示差分格式的解圖3-1、圖3-3對比可以看出，精確解和古典顯格式（網(wǎng)格比r=0.5）的圖形是一致的。圖3-4精確數(shù)值解、古典顯格式的Y-Z平面圖結果可以看出古典顯格式的邊界值和精確解一樣。圖3-5是r分別取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 古典隱格式close allclcT=0.2X=1.0M=41N=21ra=(T/(M-1)/(X/(N-1)2); %網(wǎng)格比fprintf(穩(wěn)

13、定性系數(shù) S=ra 為:n);disp(ra); %顯示網(wǎng)格比u=zeros(M,N); %構造一個M行N列的矩陣用于存放時間t和變量x for i=2:N-1 xx=(i-1)*(X/(N-1); u(1,i)=sin(pi*xx); %t=0時刻的賦值，邊界條件 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)對矩陣A進行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的上次對角線即為A的上次對

15、角線 U(i,i)=A(i,i)-L(i,i-1)*U(i-1,i); endfor k=1:M-1 %外層大循環(huán)是求時間網(wǎng)格2,3，M的求解u%以下是追趕法的求解過程%-追的過程-即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%-趕的過程-即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 %追趕法結束 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)格對每個點賦值再畫圖surf(x,t,u);xlabel(x),ylabel(t),zlabel(u);title(古典隱格式); %此程序得到圖是3-7圖3-7古典隱格式程序結果圖（r=2）圖3-8精確數(shù)值解、古典隱格式程序結果的Y-Z平面圖（r=2）圖3-9古典隱格式在取不同網(wǎng)格比時的誤差傳播結果圖圖3-10 不同時間取值時精確解、與古典隱格式的值對比圖（網(wǎng)格比r=2）

17、紅線表示精確解、藍色線表示差分格式的解圖3-7古典隱格式在r=2的圖形與精確解是一致的。圖3-8精確數(shù)值解、古典隱格式的Y-Z平面圖結果可以看出古典隱格式在t=0.2處的值的邊界值和精確解還是有誤差的。圖3-5是r分別取0.245、0.5、0.72、1.125時的誤差傳播圖像，邊界位置網(wǎng)格數(shù)為5處的誤差為0.01得到的，可以看出r取不同的值時都是穩(wěn)定的；即古典隱格式對任意的網(wǎng)格比穩(wěn)定性。從圖3-10可以看出隱格式隨著時間的增加，差分格式計算的結果和精確結果越來越大；隱格式雖然對任意網(wǎng)格比都是穩(wěn)定的，但是計算的精確度是它的不足。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); %構造一個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;%因為Richadson格式需要知道前兩行的值，第二行值我們采用隱格式求得%下面是隱格式求解第二行，和隱格式的程序一樣，只是求一行，此處不再做注釋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的上次對角線即為A的上次對角線 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%通過隱格式已求得第二行的值u(2,i)，下面是Richardson格式的過程for k=2:M-1 %時間軸第3行開始到第M行 for i=2:N-1 %位置網(wǎng)格點 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ū)域劃分進行賦值畫圖surf(x,t,u);xlabel(x),ylabel(t),zlabel(u);title(Richardson格式); %此程序得到的圖形是圖3-11圖3-11 Richardson顯格式程序結果圖（r=0.5）圖3-12精確數(shù)值解、Richardson顯格式程序結果的Y-Z平面圖（r=0.5）圖3-13 Richardson顯格式在取不同網(wǎng)格比時都不穩(wěn)定的結果圖圖3-11是 Richardson顯格式在r=0.5時的程序結果圖，可以看

22、出不穩(wěn)定。圖3-12是精確數(shù)值解、Richardson顯格式程序結果的Y-Z平面圖。從圖3-13可以看出Richardson顯格式在取不同網(wǎng)格比時都出現(xiàn)不穩(wěn)定現(xiàn)象，和理論結果相一致。所以說Richardson顯格式絕對不穩(wě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); %構造一個M行N列的矩陣用于存放時間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格式需要兩個矩陣，下面的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; 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); %下面對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;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); %首先填充右邊向量，然后進行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的上次對角線即為A的上次對角線 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)，即時間步長 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是原來的列向量矩陣乘上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的值全部求出，以下畫圖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格式程序結果圖（r=2）圖3-15精確數(shù)值解、Crank-Nicholson格式程序結果的Y-Z平面圖（r=2）圖3-16 Crank-Nicholson格式在取不同網(wǎng)格比時的誤差傳播結果圖圖3-17 不同時間取值時精確解、與C-N格式的解對比圖（網(wǎng)格比r=2）紅線表示精

28、確解、藍色線表示差分格式的解圖3-14是程序運行得到的Crank-Nicholson格式在網(wǎng)格比取r=2的結果，和精確解圖像一致。在t=0.2時從圖3-15 的Y-Z面可以看出結果還是有一定的誤差。理論上Crank-Nicholson格式對任意的網(wǎng)格比也是穩(wěn)定的，從圖3-16可以看出在r取0.245、0.5、0.72、1.125誤差傳播圖像可以看出誤差不會擴撒。圖3-17是不同時間取值時精確解、與C-N格式r=2時的解對比圖，計算還是具有誤差。 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); %構造一個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 %其實初始矩陣已經(jīng)將i=1和i=N列的初值賦零了 u(k,1)=0; u(k,N)=0; end;%第二行用隱格式求得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; %第一行第二個和最后一行倒數(shù)第二個一個賦值 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的上次對角線即為A的上次對角線 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格式); %此程序為Richardson格式的改進，得到圖3-18圖3-18 Du Fort Frankle格式程序結果圖（r=2）圖3-19精確數(shù)值解、Du Fort Frankle格式程序結果的Y-Z平面圖（r=2）圖3-20 Du Fort Frankle格式在取不同網(wǎng)格比時的誤差傳播結果圖圖3-21 不同時間取值時精確解、與D-F格式的解對比圖（網(wǎng)格比r=2）紅線表示精確解、藍色線表示差分格式的解D-F格式是Richardson格式（絕對不穩(wěn)定）的改進，從圖3-18可以看出當r時D-F格式是穩(wěn)定的；圖

33、3-20表示D-F格式網(wǎng)格比r為0.245、0.5、0.72、1.125時誤差傳播圖像，不同的網(wǎng)格比誤差都不會擴撒。說明這種格式是穩(wěn)定的，和理論上的結果相一致。圖3-21是不同時間取值時精確解、與D-F格式的解對比，隨時間的增加計算的值和差分得到的值有誤差。此種格式雖然是絕對穩(wěn)定的，但是計算的精度還是有待提升。四、結論及感想從程序得出的結果與精確解的對比來看和理論是上的結果基本一致。比如古典顯格式網(wǎng)格比r小于等于0.5（偏微分方程的系數(shù)a取1）才穩(wěn)定，從MATLAB編程運行的結果也可以看出r小于等于0.5是穩(wěn)定的，r大于0.5不穩(wěn)定。又如Richardson格式理論上是絕對不穩(wěn)定的，從編程的結

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

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

36、的對我們理解有限差分法和求解方程提供了很大的幫助。最后，感謝張老師在課堂上的知識的講解；同時也感謝在完成報告期間對我提供幫助的同學！附 錄誤差傳播三維圖（以顯格式（圖3-5）為例）：close allclcT=0.2X=1.0M=41N=8u=zeros(M,N); %構造一個M行N列的矩陣用于存放時間t和變量xra=(T/(M-1)/(X/(N-1)2); %網(wǎng)格比fprintf(穩(wěn)定性系數(shù) S=ra 為:n);disp(ra); % 顯示網(wǎng)格比u(1,5)=0.01;%即t=0時刻賦值，誤差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。由行開始 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)格對每個點賦值再畫圖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時刻賦值，誤差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時刻賦值，誤差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