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

1、南京理工大学课程考核论文课程名称：高等数值分析论文题目：有限差分法求解偏微分方程姓名：罗晨学号：1成绩：任课教师评语：签名：年 月日有限差分法求解偏微分方程、主要内容1. 有限差分法求解偏微分方程，偏微分方程如一般形式的一维抛物线型方程:2uf (x,t)（其中为常数）6 / 30具体求解的偏微分方程如下:x2u(x,O) sin( x)u(O,t) u(1,t)0 t 02. 推导五种差分格式、截断误差并分析其稳定性；3. 编写MATLAB程序实现五种差分格式对偏微分方程的求解及误差分析;4. 结论及完成本次实验报告的感想。、推导几种差分格式的过程:有限差分法（finite-differen

2、ee method是一种数值方法通过有限个微分方 程近似求导从而寻求微分方程的近似解。有限差分法的基本思想是把连续的定解 区域用有限个离散点构成的网格来代替； 把连续定解区域上的连续变量的函数用 在网格上定义的离散变量函数来近似；把原方程和定解条件中的微商用差商来近 似，积分用积分和来近似，于是原微分方程和定解条件就近似地代之以代数方程 组，即有限差分方程组，解此方程组就可以得到原问题在离散点上的近似解。推导差分方程的过程中需要用到的泰勒展开公式如下：f(x)f(Xo)(! ) (x Xo)2|)(x Xo)(o)(X Xo) o(x Xo) n!求解区域的网格划分步长参数如下：(2-1)tk

3、 1 tkxk ixkh(2-2)2.1古典显格式2.1.1古典显格式的推导由泰勒展开公式将u(x,t)对时间展开得U(冷t)u(Xi,tk)叶)i,k(ttk) o(t tk)2)(2-3)当t tk i时有u(Xi,tk i)u(X,tk)U(Xi,tk)u()i,k(tk(寸)i,ktk) o(tk1 tk)2)(2-4)0(2)得到对时间的一阶偏导数u()i,k =由泰勒展开公式将u(x,t)对位置展开得u(Xi,tk i)u(Xi,tk)0(2-5)u(X,tk)u(Xi,tk)ui()i,k(XXi)(x2!2u2)i,k(XXi)2X3o(x Xi)(2-6)当x x 1禾口 x

4、Xi 1 时,代入式(2-6)得u(X i,tk)u(Xi,tk)(f(Xi1Xi)u(X i，tju(Xi，tk)(如 XiiXi)1丄2!2u2)i,k(x 1 Xi)o(Xi 1X2口r)i,k(x 1 X)20( X 1xx)3)X)3)(2-7)因为Xk 1 Xkh，代入上式得u(n i，tk) u(A，tQH)i,kXu(N i,tk) u(Xi,tk)(丄)i,kX£(j)i,k h2 o(h3)2! x£(理)i,k h2 o(h3)2! x(2-8)得到对位置的二阶偏导数2/ uudiitk)gjjlxAu(x i，tk)小2、(2 )i,k, 2o(h

5、)Xh2(2-9)将式(2-5)、(2-9)代入一般形式的抛物线型偏微分方程得u(Xi,tk i) u(X,tQu(Xi,tk) 2u(Xi,tk) u(Xi1,tk)f(Xi,tk) o(h2h2)(2-I0)为了方便我们可以将式(2-10)写成(2-11)k 1 kUiUikk kUi 1 2UiUi 1h2k 1 kUi Ui2ukkUi 1最后得到古典显格式的差分格式为uk 1(1 2ra)ukkUi 1kUi 1其中r孑，古典显格式的差分格式的截断误差是o(h2) o2.1.2古典显格式稳定性分析古典显格式(2-13)写成矩阵形式为k 1Uh1 2ra IraCkUh(2-12)(2

6、-13)(2-14)k z k kkUh(U1 ,U2,Unk2 ,UN1)01C 0M000M10 (N 1) (N 1)上面的C矩阵的特征值是:H 1 2ra IraCH 1 2rara c =:12ra2ra cos( jh)1 2ra 1cos( jh)1 4ra sin2 j h2j1,2,N即1 14rasin 丿h 120ra 12cos( j h)C12使(H)1 ,j 1,2,(2-15)结论：当0 ra11时'所以古典显格式是稳定的。2.2古典隐格式8 / 302.2.1古典隐格式的推导将t 代入式(2-3)得u(Xj,tkJ u(Xj,tk)(十)j,k(tki

7、tk) o(tk ! tk)2)u(Xj,tki) u(Xj,tk)(汁)j,k0(2)得到对时间的一阶偏导数u u(Xj,tk) u(Xj,tk i)(-)j,k=0()将式(2-9)、(2-18)原方程得到u(Xj,tk) u(Xj,tk i)u(Xj i,tk) 2u(Xj,tk) u(Xj i,tk)为了方便把(2-i9)写成k k iUj Ujkk kUj i 2uj Uj ih2k k ikk kuj uj2 uj i 2uj uj ih最后得到古典隐格式的差分格式为c kkkk i(i 2ra)ujr uj i uj i uj其中r -，古典隐格式的差分格式的 截断误差是o( h

8、2) h2.2.2古典隐格式稳定性分析将古典隐格式(2-22)写成矩阵形式如下i 2ra I raC uhk i ujfhk(r r)h误差传播方程i 2ra I raC vj i vjA i 2ra I raC , B I所以误差方程的系数矩阵为iiH A i 2ra I raC(2-i6)(2-i7)(2-i8)o(h2)(2-i9)(2-20)(2-2i)(2-22)(2-23)(2-24)12 / 301 2ra 2racosj hj 1,2,N 1使(H)1，显然11 2ra 2 ra cos( j h)11 2ra (1 cos( j h)12 j h1 4ra sin21恒成立。

9、结论：对于r 0,即任意网格比下，古典隐格式是绝对稳定的2.3 Richardson 格式2.3.1 Richards on格式的推导将t tk 1和t tk 1，代入式(2-3)得u(N,tk1)u(Xi,tk) (U)i,k(tk 1 tk) o(tk 1 tk)2)t(2-25)u2u(x,tkju(Xi,tk) ()i,k(tk1 tk) o(tk 1 tk)即u(冷 tk 1)u(N,tk)(十)i,ko( 2)u(x,tk1)u(x,tk) (U)i,ko( 2)(2-26)由此得到可得)i,kU(X,tk 1) U(X,tk 1)2o( 2)(2-27)将式(2-9)、(2-27

10、)代入原方程得到下式U(Xi ,tk 1) u(x ,tk 1)2u(Xi 1,tk) 2u(x,tk) u(Xi1,tk)h2f(X,tk) o( 2 h2)(2-28)(2-29)为了方便可以把式(2-28)写成k 1 k 1UiUi2kk kUi 1 2Ui Ui 1k 1Uik 1Ui最后得至U Richards onkk k-r Ui 1 2q Ui h显格式的差分格式为(2-30)k 1Uik k k2r Ui 1 2Ui Ui 1k 1Ui2 fk(2-31)其中r -，古典显格式的差分格式的 截断误差是o( 2 h2)h2.3.2 Richards on稳定性分析将Richar

11、dson显格式(2-31)写成如下矩阵形式1uk 12r C2IkUh2 fhk(2-32)误差传播方程矩阵形式k 1VhkVh2rkVh2IkVhk 1Vh(2-33)再将上面的方程组写成矩阵形式k 1Vk-V2ra(CI2I) I0(2-34)系数矩阵的特征值是4ra cos(j1h)4ra解得特征值为1,2Max 18ra sin2(2-35)2 j h22.4 j h8ra sin . -64r a sin 42 Y222 =4ra sin2+J+16r2a2sin42 21 (恒成立)(2-36)(2-37)结论：上式对任意的网比都恒成立，即Richardson格式是绝对不稳定的13

12、 / 304. Crank-Nicholson 格式3.4.1 Crank-Nicholson 格式的推导 将t tk i代入式(2-9)得u(Xj,tk 1) u(Xj,tk) (-U)j,k(tk； tk) o(tk 1 tk)2)t(2-40)u2u(Xj,tk 2) u(Xj,tk+i) ()j,ki(tkg tk+i) o(tk 1 tk+i)即u2u(Xj,tk i) u(Xj,tk) ( -)j,k - o()t 2(2-4i)u(Xj,tk i) U(Xj,tk+i)(汁)j,k i - o( 2)得到如下方程u 2 u(Xj,tk i) u(Xj,tk)2(下爪二-o()(2

13、-42)u2u(Xj,tki)u(Xj,tki)2(-)j,ki=-o()所以t tk i处的一阶偏导数可以用下式表示：N 2i uui 2 u(Xj,tk i) u(Xj,tk)2 u(Xj,tk i) u(Xj,tk i)2-()j,k i ( -)j,k-o()u(Xj,tk i) u(Xj,tk)“ 2、o()(2-43)将t tk ,X Xj i和x Xj i代入式(2-6)可以得到式(2-9)；同理t tk i ,x Xj i和x Xj i代入式（2-6）可以得到2(2-44)uu(Xj i,tk i) 2u(Xj,tk i) u(Xj i,tk i) 亠2、()j,k i2o(h

14、 )xh所以t tk i，Xj处的二阶偏导数用式（2-6）、（2-44）表示:1 u(Xj1,tk1)2u(Xj,tk1)u(Xj1,tk1)u(Xj 1,tk)2u (Xj,tk)u (Xj 1,tk)o(h2)2h2h2(2-45)所以t tk 1，Xj处的函数值可用下式表示：12 f (Xj,tk 1)f (Xj,tk)(2-46)原方程变为：1 ( U)2 ( t)j,k12 2 1(*)j,k(2 )j,k (2)j,k 1-t2 xx2fjk1 f?o(h2)(2-47)2(上)(_u)2 j,k2 j,k 1xx将差分格式代入上式得:u(Xj,tk 1) U(Xj,tk)U(Xj

15、1,tk 1) 2U(Xj,tk1) U(Xj 1,tk 1)U(Xj 1,tk)2U(Xj,tk) U(Xj 1,tk)2h2h21f (Xj,tk 1)22f(Xj,tk)o(h2)为了方便写成：(2-48)k 1k rk 1UjUj uj !k 12Ujk 1kkUj 1Uj 1 2Ujk1_ k 1U 12 fjfjk(2-49)2 2最后得到Crank-Nicholson格式的差分格式为/k 1 r k 1 k1 k r k k1- k 1- k tr.1 r UjUj 1 Uj 1 = 1 r Uj Uj 1 Uj 1 +fjfj (2-50)其中r , Crank-Nichols

16、on格式的差分格式的 截断误差是o( 2 h2)。h341 Crank-Nicholson 稳定性分析Crank-Nicholson格式写成矩阵形式如下:rk 1(1 r )I 2 C 山=rk(1 r )I 2 C Uh +1 k 1k2 fhfh(2-51)误差传播方程是：(1 r )I；Ck 1Vh = (1 r )IrkC Vh2(2-52)可以得到：(1 r )1(1 r )IH (1 r ) ra cos(j h)(1 r ) ra cos(j h)2 j h1 2ra sin21 2rasin2jJh2使(H)1即1 2ra sin2-1 九11 2ra sin2 -22 j h

17、2 j h2 j h1 2ra sin1 2ra sin1 2ra sin2 2 22 j h2 j h1 2ra s in1 2ra s in2 22 j h2 j h1 2ra sin1 2ra sin2 22rasin2U 2rasin2U2 2rasing2rasin22(2-53)(2-54)(2-55)(2-56)(2-57)(2-58)上式恒成立。结论：Crank-Nicholson格式对任意网格比也是绝对稳定的5. Du Fort Fran kle 格式(Richards on 格式的改进)1将Uik -(Uik 1 uk 1)代入式(2-31)并化简得到 Du Fort Fr

18、ankle：uk 1 2ru'1uk 1uk 1u'1 uk 1 2(2-59)(1 2ra)uk12ruik1uik1(12ra)u： 12 fik(2-60)可以证明Du Fort Frankle格式是绝对稳定的。因为此格式是Richardson格式的改 进格式，因此截断误差还是o( 2 h2)。3.6总结(1) 古典显格式古典显格式的差分格式为uk 1 (1 2ra)u： r u uk1fik截断误差：o(h2)。1稳定性：当0 ra -时，古典显格式稳定。2(2) 古典隐格式古典隐格式的差分格式为kkkk 1k(1 2ra)Uj r Uj 1 Uj 1Ujf j截断误差

19、：o(h2) o稳定性：任意网格比古典隐格式绝对稳定。Richards on显格式Richardso n显格式的差分格式为k 1ck /* k kk 1r kui2r ui 1 2ui ui 1 u 2 fi截断误差：o( 2 h2) 稳定性：任意网格比Richardson格式绝对不稳定 Crank-Nicholson 格式Crank-Nicholson格式的差分格式为/k 1 r k 1 k1 k r k k1- k 1- k1 r Uj y Uj 1 Uj 1 = 1 r Uj y uj 1 uj 1 + - fjfj截断误差：o( 2 h2) o稳定性：Cran k-Nicholson格

20、式对任意网格比绝对稳定。 Du Fort Frankle 格式(1 2ra)u：12ru uik1(1 2ra)uk 1 2 fik 截断误差：o( 2 h2) o稳定性：Du Fort Fran kle格式对任意网格比绝对稳定三、MATLAB实现五种差分格式对偏微分方程的求解及误差分析3.1精确数值解u(x,O) sin( x)u(O,t) u(1,t)0t 0上述偏微分方程的精确解是2tu(x,t) e sin( x) (0 x 1, t 0) 区域取值范围：0 x 1 , 0 t 0.2。用MATLAB对精确解进行编程画出三维图像精确解程序如下:close allclcx,t=meshg

21、rid(0:0.01:1,0:0.001:0.2)u=exp(-pi*pi*t).*si n( pi*x)mesh(x,t,u);surf(x,t,u);xlabel('x'),ylabel('t'),zlabel('u');title('精确数值解');shadi ng in terp0.80.6L.0502图3-1精确数值解的三维图24 / 30(a)精确数值解X-Y平面(b)精确数值解X-Z平面(c)精确数值解Y-Z平面图3-2精确数值解的三个平面图3.2五种差分格式MATLAB程序3.2.1古典显格式close allcl

22、cT=0.2X=1.0M=41N=11u=zeros(M,N);%构造一个M行N列的矩阵用于存放时间t和变量xra=(T/(M-1)/(X/(N-1)A2);% 网格比fprintf('稳定性系数S=ra 为:n');disp(ra);%显示网格比for i=2:N-1xx=(i-1)*(X/(N-1);u(1,i)=si n( pi*xx);end;%即t=0时刻赋值，边界条件for k=1:Mu(k,1)=0;u(k,N)=0;en d;% x=0,x=1处的边界条件for k=1:M-1%矩阵是从y轴表示行k，x轴表示列i。由行开始for i=2:N-1u(k+1,i)=

23、(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);%将区域划分成网格对每个点赋值再画图surf(x,t,u);xlabel('x'),ylabel('t'),zlabel('u');title('古典显格式');%此程序得到的是图3-30 0图3-3古典显格式程序结果图（r=0.5）图3-4精确数值解、古典显格式程序结果的Y-Z平面图（r=0.5）26 / 30古咄显格式，逞老椿播1=0.2

24、45ii.t y r W0.D'DOS 6D古即显帝貳i£差特J& i=0Z2图3-5古典显格式在取不同网格比时的误差传播结果图34 / 30图3-6不同时间取值时精确解、与古典显格式的值对比图（网格比r=0.5）红线表示精确解、蓝色线表示差分格式的解图3-1、图3-3对比可以看出，精确解和古典显格式（网格比r=0.5）的图形是一致的。图3-4精确数值解、古典显格式的Y-Z平面图结果可以看出古典显格 式的边界值和精确解一样。图 3-5是r分别取0.245、0.5、0.72、1.125时的误差 传播图像，边界位置网格数为 5处的误差为0.01得到的，可以看出r小于等于

25、0.5是稳定的；但是r大于0.5出现不稳定现象。很好的验证了古典显格式稳定性。3.2.2 古典隐格式close allclcT=0.2X=1.0M=41N=21ra=(T/(M-1)/(X/(N-1)A2);fprintf(' 稳定性系数 S=ra disp(ra);u=zeros(M,N);for i=2:N-1 xx=(i-1)*(X/(N-1); u(1,i)=sin(pi*xx);end;for k=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;%网格比 为 :n'

26、);%显示网格比%构造一个 M 行 N 列的矩阵用于存放时间%t=0 时刻的赋值，边界条件% x=0,x=1 处的边界条件 %隐格式的矩阵形式中的 A 矩阵赋值t 和变量 xA(1,2)=-ra;A(N-1,N-2)=-ra; for m=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,1);for i=1:N-1d(i,1)=sin(pi*(i-1)*(

27、1.0/(N-1);end%隐格式右边矩阵赋值%以下循环对矩阵 A 进行 LU 分解U(1,1)=A(1,1);for i=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);endfor k=1:M-1%外层大循环是求时间网格2,3,M的求解u%以下是追赶法的求解过程% 追的过程 即 Ly=d 的求解 yy(1,1)=d(1,1);for i=2:ny(i,1)=d(i,1)-L(i,i-1)*y(i-1,1);end% 赶的过程 即 Ux

28、=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: nu(k+1,i)=x(i,1)endd=zeros(N-1,1);%更新右边矩阵d=x%每次外循环更换右边矩阵endfor k=1:Mu(k,1)=0;en d;disp(u);%显示差分法求得的值x,t=meshgrid(1:N,1:M);%将区域划分成网格对每个点赋值再画图surf(x,t,u);xlabel('x'),ylabel('t'),zlabel

29、('u');title('古典隐格式');%此程序得到图是3-7图3-7古典隐格式程序结果图（r=2）lUEfTlK* II- V.叶格it图3-8精确数值解、古典隐格式程序结果的Y-Z平面图（r=2）图3-9古典隐格式在取不同网格比时的误差传播结果图w.牯巴馭、占营-，琴書曲鹦刃比/0 15ow/吗7V&臥搭寸泊銘fff!忙Thkr=2）红线表示精确解、图3-10不同时间取值时精确解、与古典隐格式的值对比图（网格比 蓝色线表示差分格式的解图3-7古典隐格式在r=2的图形与精确解是一致的。图3-8精确数值解、古 典隐格式的丫-Z平面图结果可以看出古典隐格

30、式在t=0.2处的值的边界值和精确 解还是有误差的。图3-5是r分别取0.245、0.5、0.72、1.125时的误差传播图像,边界位置网格数为5处的误差为0.01得到的，可以看出r取不同的值时都是稳定的；即古典隐格式对任意的网格比稳定性。从图3-10可以看出隐格式随着时间 的增加，差分格式计算的结果和精确结果越来越大； 隐格式虽然对任意网格比都 是稳定的，但是计算的精确度是它的不足。3.2.3 Richardson 显格式close all clcT=0.2X=1.0000M=41N=11 ra=(T/(M-1)/(X/(N-1)A2);fprintf(' 稳定性系数 S=ra 为

31、:n');disp(ra);%构造一个 M 行 N 列的矩阵u=zeros(M,N);for i=2:N-1xx=(i-1)*(X/(N-1);u(1,i)=sin(pi*xx);%边界条件end;for k=1:Mu(k,1)=0;u(k,N)=0;% 边界条件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;for

32、m=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);for i=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);endy(1,1)=u(1,1);for i=2:ny(i,1)=u(1,i)-L(i,i-1)*y(i-1,1);endx(n,1)=y(n,1)/

33、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);endu(2,1)=0;for i=2:N-1u(2,i)=x(i,1)end%通过隐格式已求得第二行的值u(2,i)，下面是Richards on格式的过程for k=2:M-1for i=2:N-1%时间轴第3行开始到第M行%位置网格点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);%区域划分进行赋值画

34、图surf(x,t,u);xlabel('x'),ylabel('t'),zlabel('u');title('Richardson格式');%此程序得到的图形是图3-11图3-11 Richards on显格式程序结果图(r=0.5)D.D5.!36 / 30图3-12精确数值解、Richardson显格式程序结果的 Y-Z平面图(r=0.5)图3-13 Richards。n显格式在取不同网格比时都不稳定的结果图图3-11是Richardson显格式在r=0.5时的程序结果图，可以看出不稳定。 图3-12是精确数值解、Richa

35、rdson显格式程序结果的Y-Z平面图。从图3-13可 以看出Richardson显格式在取不同网格比时都出现不稳定现象，和理论结果相一 致。所以说Richardson显格式绝对不稳定，这种差分格式不能使用。后面有改进 的格式D-F格式。3.2.4 Crank-Nicholson 格式close allclcT=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);for i=2:N-

36、1xx=(i-1)*(X/(N-1);u(1,i)=sin(pi*xx);% 边界条件en d;for k=1:Mu(k,1)=0;u(k,N)=0;%边界条件en d;%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-2A(m,m-1)=-ra/2;A(m,m)=1+ra;A(m,m+1)=-ra/2;end;B=zeros(N-1,N-

37、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-2B(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:nL(i,i-1)=A(i,i-1)/U(i-1,i-1);U(i-1,i)=A(

38、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%按层外围大循环，即时间步长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(i,i)*d(i,1)+B(i,i+1)*d(i+1,1);

39、end%以上是右边矩阵的填充更新%追求解Ly=b，其中b是原来的列向量矩阵乘上B系数矩阵得到的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:1x(i,1)=(y(i,1)-U(i,i+1)*x(i+1,1)/U(i,i);endfor i=1:nu(k+1,i)=x(i,1)endd=zeros(N-1,1);%右边向量d=xendfor k=1:M u(k,1)=0;en d;disp(u); % u的值全部求出，以下画图x,t=m

40、eshgrid(1:N,1:M);surf(x,t,u);xlabel('x'),ylabel('t'),zlabel('u');title(' Crank-Nicholson 格式');%此程序得到的图像是图3-14CrankNichDlhinti1 衣图3-14 Crank-Nicholson格式程序结果图（r=2）图3-15精确数值解、Crank-Nicholson格式程序结果的 Y-Z平面图（r=2）42 / 30"ra-1 Ni-nnlL .mtt 式 IS jB rD 坤图3-16 Crank-Nichols

41、on格式在取不同网格比时的误差传播结果图I-：, m-I 崎蹄 Cm Ft Hl' U!'.g? IH g D8ei琦討即堵曜解-aittffiliwtt 1'>-£第94XDB图3-17不同时间取值时精确解、与C-N格式的解对比图（网格比 r=2）红线表示精确解、蓝色线表示差分格式的解图3-14是程序运行得到的Crank-Nicholson格式在网格比取r=2的结果，和 精确解图像一致。在t=0.2时从图3-15的丫-Z面可以看出结果还是有一定的误 差。理论上Crank-Nicholson格式对任意的网格比也是稳定的，从图 3-16可以看 出在r取0.

42、245、0.5、0.72、1.125误差传播图像可以看出误差不会扩撒。图3-17是不同时间取值时精确解、与 C-N格式r=2时的解对比图，计算还是具有误差。3.2.5 Du Fort Frankle 格式 close all clcT=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 列的矩阵for i=2:N-1xx=(i-1)*(X/(N-1);u(1,i)=sin(pi*xx);%i 表示 x 轴，k 表示 y

43、轴(即 t)end;for k=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)=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-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,

44、1); U(1,1)=A(1,1);for i=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% 追 y(1,1)=u(1,1);for i=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);for i=n-1:-1:1x(i,1)=(y(i,1)-U(i,i+1)*x(i+1,1)/U(i,i);endu(2,1)=0;for i=2:N-1u(2,i)=x(i,1

45、)end%第二行求出，下面用D-F 差分格式for k=2:M-1for i=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(1:N,1:M);surf(x,t,u);xlabel('x'),ylabel('t'),zlabel('u');3-18title('Du Fort Frankle格式');%此程序为 Richardson格式的改进，得到图Dv Fon式0 口图3-1

46、8 Du Fort Frankle格式程序结果图（r=2）1E确剤值解D Fort F (flnk fir FC43 / 30图3-19精确数值解、Du Fort Frankle格式程序结果的 Y-Z平面图（r=2）图3-20 Du Fort Frankle格式在取不同网格比时的误差传播结果图图3-21不同时间取值时精确解、与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 时误差传播图像，不

47、同的网格比误差都不会扩撒。说明这种格式是稳定的，和理论上的结果相一致。图3-21是不同时间取值时精确解、与 D-F格式的解对比， 随时间的增加计算的值和差分得到的值有误差。此种格式虽然是绝对稳定的，但是计算的精度还是有待提升。四、结论及感想从程序得出的结果与精确解的对比来看和理论是上的结果基本一致。 比如古 典显格式网格比r小于等于0.5（偏微分方程的系数a取1）才稳定，从MATLAB 编程运行的结果也可以看出 r 小于等于 0.5 是稳定的， r 大于 0.5 不稳定。又如 Richardson 格式理论上是绝对不稳定的，从编程的结果在取不同的网格比 Richardon 格式都是不稳定的，理

48、论和结果一致。理论上对不稳定格式进行改进 使其稳定，比如得到的 D-F 格式，取不同的格式网格比都是稳定的，很好的验 证了改进的格式的稳定性。本次报告首先推导了一维抛物线型偏微分方程的一般差分格式。 分别是古典 显格式、古典隐格式、Richards on显格式、C-N格式进版的Richards on格式D-F 格式。推导中得到各种格式的截断误差， 从理论上分析了各种格式的稳定性。 对 于不稳定的格式进行修改得到稳定的格式，即 Richardson格式的修改。通过 MATLAB 编程实现了各种格式的程序实现。用实验的方法来验证理论 结果，能更好的理解各种差分格式的稳定性及操作过程。 报告中的程序

49、都是基本 程序，误差图与二维 x-u 的图像（见附录）都是在基本程序的基础上对参数的修 改得到的图像。通过这次报告的完成学到了很多的东西。首先，对编程有了进一步的了解； 尤其是使用 MATLAB 编程。在这个过程中也遇到了很多的问题，通过查阅资料 并利用网络资源寻求解决办法。 其次，对于差分格式在程序上的实现并不是简单 的书写，需要步步衔接好； 比如好几种格式都用到差分格式的矩阵形式， 尤其是 隐式格式不能直接求解， 需要应用追赶法进行求解； 编写过程中通过直接求解得 出的结果都不正确， 通过追赶法才能得到正确的结果。 最后， 我们专业是电磁场 与微波技术且偏计算， 经常遇到偏微分方程的求解；

50、 所以这次试验毫无疑问的对 我们理解有限差分法和求解方程提供了很大的帮助。最后，感谢张老师在课堂上的知识的讲解； 同时也感谢在完成报告期间对我 提供帮助的同学！47 / 30附录3-5）为例）：误差传播三维图(以显格式(图 close all clcT=0.2X=1.0M=41N=8u=zeros(M,N);%构造一个 M 行 N 列的矩阵用于存放时间 t 和变量 xra=(T/(M-1)/(X/(N-1)A2);% 网格比fprintf(' 稳定性系数 S=ra 为 :n'); disp(ra); % 显示网格比 u(1,5)=0.01;%即 t=0 时刻赋值，误差0.01f

51、or 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-1u(k+1,i)=(1-2*ra)*u(k,i)+ra*u(k,i+1)+ra*u(k, i-1);%此处为古典显格式endend disp(u); % 显示差分法求得的值 x,t=meshgrid(1:N,1:M);%将区域划分成网格对每个点赋值再画图 subplot(2,2,1) surf(x,t,u);xlabel('x'),ylabel('t'),zlabel('u');title(' 古典显格式 -误差传播 r=0.245'); clcT=0.2X=1.0M=41N=11 u=zeros(M,N); ra=(T/(M-1)/(X/(N-1)A2);fprintf(' 稳定性系数 S=ra 为 :n'