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

1、南京理工大学课程考核论文课程名称：高等数值分析论文题目:有限差分法求解偏微分方程姓 名：罗 晨学 号： 115104000545成 绩：任课教师评语：签名：年月日有限差分法求解偏微分方程一、主要内容1. 有限差分法求解偏微分方程，偏微分方程如一般形式的一维抛物线型方程:2 = f (x,t)(其中 a为常数).:t .x具体求解的偏微分方程如下:_2Q<x<1；u : uf 2 = 0.t .xu(x,0) =sin(二 x)u(0,t) =u(1,t) =0t _0(2-2)2.1古典显格式2. 推导五种差分格式、截断误差并分析其稳定性；3. 编写MATLAB程序实现五种差分格式

2、对偏微分方程的求解及误差分析;4. 结论及完成本次实验报告的感想。二、推导几种差分格式的过程:有限差分法（finite-difference methods是一种数值方法通过有限个微分方 程近似求导从而寻求微分方程的近似解。有限差分法的基本思想是把连续的定解 区域用有限个离散点构成的网格来代替； 把连续定解区域上的连续变量的函数用 在网格上定义的离散变量函数来近似；把原方程和定解条件中的微商用差商来近 似，积分用积分和来近似，丁是原微分方程和定解条件就近似地代之以代数方程 组，即有限差分方程组，解此方程组就可以得到原问题在离散点上的近似解。(n) ,一(x x0)n o(x x°)n

3、 n!(2-1)推导差分方程的过程中需要用到的泰勒展开公式如下：f（x） = f （x0）（乂一乂0） f 2：" （x -x。）2求解区域的网格划分步长参数如下：tk 1 -tk = xk 1 F =h2.1.1古典显格式的推导由泰勒展开公式将u(x,t)对时间展开得u(x,tX u(iXktw£(,i k 项 2)。顶kt ) )(2-3)；:t当t =tk中时有U，2u( x, ti X u( x,kt ) ( , A. tktk ) oi(tk )(2-4)=u( x, k )半，i)k , o2()ft得到对时间的一阶偏导数(当尸心心一心小项,(2-5);t由泰勒

4、展开公式将u(x,t)对位置展开得-U、，、1 , f 2U、，、2、3、,c C、u(x,tk) =u(xi,tk) ()i,k(x-xi)(r)i,k(x xi)o(x-x) )(2-6)x2! :x当x =为中和x =为时，代入式(2-6)得u1 ,；：2u23、u(x+,tk) =u(xi,tk)十(丁)i,k(x+xi)十指(>)i,k(x4i xi) +o(x xi)-2(2-7):x2!:x；:u1 ；：2u23u(",tk) =u(xi,tk)+(=)i,k(x口一) m)i,k(x_ixi) +o(xj x)2! : x因为xz -xk = h ,代入上式得u

5、、u(x i,tk) =u(x,tQ ( )i,k :xu、 u(Xjtk) =u(x,tk) -(一)筋 x-221- u 23h(2)i,k h2 o(h3)2! :x 2 一 一h ；(-y)i,k h o(h )2! :x(2-8)2.2古典隐格式得到对位置的二阶偏导数-2(2-9)h2(徐)_u(xf,tk) 2u(xi,tk) + u(x,tk) + o(h2)：'x ,将式(2-5)、(2-9)代入一般形式的抛物线型偏微分方程得u(xi,tk i)u(x,tk) 一 u(x i,tk)2u(xi,tk) u(xi,tk)2.- U = f(xi,tk) o( - h )-

6、h为了方便我们可以将式(2-10)写成k 1 kkk k(2-11)(2-12)Ui -UiUii-2Ui Ui 1 k-2= fi一 hk 1 k ' Ji. kk kkUi一 Ui 2 Ui 1-'2ui U=.fih最后得到古典显格式的差分格式为u：* = (12ra)Uik +ra (u% +U，)+t f：(2-13)其中r=：，古典显格式的差分格式的 截断误差是o(T+h2)。h2.1.2古典显格式稳定性分析古典显格式(2-13)写成矩阵形式为(2-14)U：1 = 1 -2ra I raC u：fhk其中 r , u： =(uk,u：,.,uN/,uN_l)。h一

7、0 1 III01 0 10C = 0：:1 0 1-0川 1。一 (N 欢 NF上面的C矩阵的特征值是：农=2cos(怦h)j =1,2,., N -1H = 1 -2ra I raCb = 1 - g ra，C = 1a?ra2 cjo S ()= 12ra - cojS( )(2-15)=1 - 4a si2nj =1, 2 ,N = . . . ,12使 P(H ) M1 ,即-1 1 - 4rasin2jlh 12。三ra -】21结论：当。壬ra1时，所以古典显格式是稳定的。22.2.1古典隐格式的推导将t=tk代入式(2-3)得u(Xj,tk)=u(Xj,tQ (兰)j*(tkf

8、) o(tkj -tk)2) ；tu(Xj,tkD=u(Xj,tk)-(m)j,ko( 2);t得到对时间的一阶偏导数(2-16)(2-17)(Cu= u(Xj ,tk) u(Xj,tk i) +ft(2-18)将式(2-9)、(2-18)原方程得到u(Xj,tk) -u(Xj,tk)u(Xj i,tk) -2u(Xj,tk) u(Xj,tk)-a ah2一， 、,一2、=f (Xj,tk) o(h )为了方便把(2-19)写成k k JUj - UjOlTu-2u：+u"l k 2 = f i _ h2 J jk k,. kk kkUj-Uj叩3fj最后得到古典隐格式的差分格式为(

9、1 2ra)u： -r： u" u' =u："(2-19)(2-20)(2-21)(2-22)其中r=：，古典隐格式的差分格式的 截断误差是o(i+h2)h2.2.2古典隐格式稳定性分析将古典隐格式(2-22)写成矩阵形式如下_ 1 2ra I -raC uk 1 = u： fk(=#)(2-23)误差传播方程1 2ra I -raC v： 1 = v：(2-24)A = 1 2ra I - raC , B = I所以误差方程的系数矩阵为.1_ =H = A - |J1 - 2ra I -raC使P(H) <1,显然11 2ra i2ra cos j 二 hj

10、 =1,2,., N -11 2ra - 2ra cos( j 二 h)11 2ra(1 -cos(j 二 h)2 j 二h1 4ra sin 2结论：对丁 Vr0,即任意网格比下，古典隐格式是绝对稳定的2.3 Richardson 格式2.3.1 Richardson格式的推导将t =tk噂和t=tk，代入式(2-3)得U(Xi,tk1) =u(xi,tk)(弋)i,k(tk 1 -tk) o(tk1 -tk)2)：t(2-25)u(Xi,tkj) =u(Xi,tk)(史)i,k(tk-tk) o(tkj -tk)2):t即f u2u(Xi,tk 1) =u(Xi,tk) ()i,k o()

11、H(2-26)；:u2u(Xi,tk)=u(x,tk) -(7)i,k , o():'t由此得到可得(2-27)（色）=u（X,tk 1）u（xi,tk Q . o（ 2）；:t i,k2将式（2-9）、（2-27）代入原方程得到下式u(xiJk Q -u(xiJk 如 一 u(xi iJD -MWJQ22、2f (XI ,tk) o( h )-h(2-28)h2最后得到Richardson显格式的差分格式为(2-29)(2-30)(2-31)(2-32)(2-33)(2-34)(2-35)(2-36)(2-37)为了方便可以把式(2-28)写成k1 kAkk kUi -UiU| -2

12、U| +Uj_a 12k 1 k d kk kkUi -Ui乏 u| i -2山 u|=.f|hk -1kk kkku|2r* i q 4 _2wu|u|2 f|其中=：，古典显格式的差分格式的 截断误差是o(/+h2) h2.3.2 Richardson稳定性分析将Richardson显格式(2-31)写成如下矩阵形式u； 1 =2r： C -2I u； u" 2 fhk误差传播方程矩阵形式v； 1 =2rC -2I W k kVh = Vh再将上面的方程组写成矩阵形式wk+='L2ra(C2I)I "kWh L苛厂J Iow系数矩阵的特征值是4ra cos(：

13、h) -4ra 1ij IL 1o22二 h 8ra sin1=02解得特征值为-8ra sin2± 64 r2a2 sin4 - + 422十2t r, 112_i，h1'2 24 j h/ L一, 、. 、 Max 仁耳， =4ra sin +1+16r a sin>1 (何成立)结论：上式对任意的网比都包成立，即 Richardson格式是绝对不稳定的4. Crank-Nicholson格式3.4.1 Crank-Nicholson 格式的推导将t =队代入式（2-9）得fm/Wz 0(")2):u2u(Xj ,tM2u(Xj, tk+1 ) *(=山-

14、书(tk 书tk+i ) *o(tk 书一tk+i)；tu(Xj,tk ,Hu(Xj,tk) Fu)j,0( 2)一 t 2u, 2、u(Xj,tk 2) =u(Xj ,tk+1) -()j,k 1 二 o()t 2(2-40)(2-41);u()j,k =:t2(u(Xj,tkG -u(Xj,tk)To( 2)(2-42)fu(T)j,k 1:t2(u(Xj,tkH2) u(Xj,tkG )0( 2)所以t =t3处的一阶偏导数可以用下式表示:1 q、.如c (一')j,k 1( ., )j,k2 I t :t2&(为儿|) u(Xj,tk)2 u(Xj,tk *) -u(X

15、j,tk.i)o( 2)u(Xj,tk i) -u(Xj,tk), 2、0()(2-43)将t =tk ,x = Xj4i和x =Xj代入式（2-6）可以得到式（2-9）;同理t =tk + ,X =Xj虫和x =Xj代入式（2-6）可以得到；：2Uu(Xj 1,tk 12u(Xj,tk 1) u(Xj J,tk 1) 儿2、(r)j,k 1 =乙o(h ):Xh2(2-44)所以t=tk , Xj处的二阶偏导数用式（2-6）、（2-44）表示:得到如下方程1 u(Xj i,tk 1)-2u(Xj,tk .1) u(Xj】,tk 1) u(Xj i,tk) -2u(Xj,tk) u(Xjj,t

16、k)=_ 2 22 ILh2h22o(h2)(2-45)所以t =tk专,Xj处的函数值可用下式表示：12 |f ( X , ' ) f ( X ,kt )(2-46)原方程变为:"(鸟 j,k (鸟 j,k12 x :x将差分格式代入上式得：u(Xj,tk 1) -u(Xj,tk):- u(Xj 1,tk 1) -2u(Xj,tk 1) u(Xj,tk i) u(Xj 1,2 2 ILh=：f (Xj,tk 书)+ f(Xj ,tk)+o(T2 +h2)1 /；：u u2 _"1 "1k .1k2= Z(fj +fj )+o(h) (2-47)tk)-2

17、u(Xj,tk) u(Xjtk)h2(2-48)为了方便写成:k 1 k rQ" k 1 k：M1 k -Uj 一山一2|luj 1 -2山jkk kuj . 12uj uj_(2-49)最后得到Crank-Nicholson格式的差分格式为< 一 k 1 r 上 k 1 k 1< 一 k r 上 k k1+. - fjk 1f jk(2-50)1 " ujuj 1 uj. = 1 T- uj ; uj 1 uj,：2其中=户，Crank-Nicholson格式的差分格式的 截断误差是o(T2+h2)。3.4.1 Crank-Nicholson 稳定性分析Cra

18、nk-Nicholson格式写成矩阵形式如下:(2-51)r一 k 1(1 r： )I 里 C uk 1= (1-r： )I2误差传播方程是:可以得到:(2-52)使P(H ) M1即H = (1+ra )眺C l' ' JTa I 今目C IC 2 J 32 )H (1- r ： ) r a co sj( h )(1 r ：卜 r a c o sj( h )1一 2ra s2口212a s i2Fjh2-1 <2 jh1 - 2a s i 目22 j 二 h1 2ra sin2<12 j： h2 j h2 j . h一 1 -2 a si _ 1 r2i sr-

19、s 1r a2s i n2222 j 二 h2 j 二 h1 -2ra sin 1 2ra sin -222 j 二 h2 j 二 h一1 2rasin -1 -2ra sin 222 j-h2 j h-2ra sin 2ra sin 222 j 二 h2 j 二 h-ra sin _ 1 一 ra sin 22(2-53)(2-54)(2-55)(2-56)(2-57)(2-58)上式包成立。结论：Crank-Nicholson格式对任意网格比也是绝对稳定的5. Du Fort Frankle 格式(Richardson格式的改进)将u =1(uik+uik_l)代入式(2-31)并化简得到

20、 Du Fort Frankle：(2-59)(2-60)k 1k k 1 k-1kk- 1 kui =2r ui 1 w -uiui 4 ui2 fi(1 2ra)uk 1 =2r： u" u， (V2ra)u 2 fik可以证明Du Fort Frankle格式是绝对稳定的。因为此格式是Richardson格式的改 进格式，因此截断误差还是oR2+h2）。3.6总结(1) 古典显格式古典显格式的差分格式为u：* =(12ra)Uik +ra (u» +心)+"，截断误差：o(i + h2) o1稳正性：当Ora壬时，古典显格式稳定。2(2) 古典隐格式古典隐格

21、式的差分格式为 kk kkk(1 2ra)Uj - r：z 凹 1 4=4. f j截断误差：o(i + h2) o稳定性：任意网格比古典隐格式绝对稳定。(3) Richardson 显格式Richardson显格式的差分格式为u：* =2ra (uL 2uk +u，)+qk，+ 2w fik截断误差：o(T2 + h2)。稳定性：任意网格比Richardson格式绝对不稳定。(4) Crank-Nicholson 格式Crank-Nicholson格式的差分格式为1 r Uk 1Uk ； uk=1 - r： Uk U " U + 1 fk 1 fk1U j2 U j 1 U jJ1

22、 .j 2 U j 1 U j j 2 J J截断误差：o( 2 h2) o稳定性：Crank-Nicholson格式对任意网格比绝对稳定。(5) Du Fort Frankle 格式(1 +2ra)Uik* =2颅(q、十口七)十(1-2佰川尸 + 27 fik截断误差：o(2 h2)0稳定性：Du Fort Frankle格式对任意网格比绝对稳定。三、MATLAB实现五种差分格式对偏微分方程的求解及误差分析3.1精确数值解一 .2；u ; u.:tu(x,0) =sin(二 x)u(0,t) =u(1,t) =0t _0图3-1精确数值解的三维图上述偏微分方程的精确解是u(x,t) =er

23、 sin(-： x)(0 £ x 壬 1, t _ 0)xlabel('x'),ylabel('t'),zlabel('u'); title('精祐敖值解');区域取值范围：0<x1 , 0壬t0.2。用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);shading interpxlabel('(a)精确数值

24、解X-Y平面(b)精确数值解X-Z平面-0.5(c)精确数值解Y-Z平面图3-2精确数值解的三个平面图3.2五种差分格式MATLAB程序3.2.1古典显格式close all clc T=0.2 X=1.0 M=41 N=11 u=zeros(M,N);%构造一个M行N列的矩阵用于存放时间t和变量xra=(T/(M-1)/(X/(N-1)A2);% 网格比fprintf('稳定仕系数 S=ra 为:n'); disp(ra);%显示网格比for i=2:N-1 xx=(i-1)*(X/(N-1); u(1,i)=sin(pi*xx); end;%即t=0时刻赋值，边界条件for

25、 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);% 此处为古典显格式end end disp(u);%显示差分法求得的值x,t=meshgrid(1:N,1:M);%将区域划分成网格对每个点赋值再画图surf(x,t,u); xlabel('x'),ylabel('t'),zlabel('u'); title(

26、9;古典显格式');%此程序得到的是图3-3图3-3古典显格式程序结果图(r=0.5)图3-4精确数值解、古典显格式程序结果的Y-Z平面图(r=0.5)图3-5古典显格式在取不同网格比时的误差传播结果图。.昭时情礁解、古成显格式的薛曲琴对比t=o 125时，情网解.古载拄格式的解曲线对比0.3£i03D.2D.40.C OS 1图3-6不同时间取值时精确解、与古典显格式的值对比图（网格比r=0.5）红线表示精确解、蓝色线表示差分格式的解图3-1、图3-3对比可以看出，精确解和古典显格式（网格比 r=0.5）的图形 是一致的。图3-4精确数值解、古典显格式的Y-Z平面图结果可以

27、看出古典显格 式的边界值和精确解一样。图 3-5是r分别取0.245、0.5、0.72、1.125时的误差 传播图像，边界位置网格数为 5处的误差为0.01得到的，可以看出r小丁等丁 0.5是稳定的；但是r大丁 0.5出现不稳定现象。很好的验证了古典显格式稳定性。3.2.2古典隐格式close allclcT=0.2X=1.0M=41ra=(T/(M-1)/(X/(N-1)A2);%网格比fprintf('稳定仕系数 S=ra disp(ra);为:n');%显示网格比u=zeros(M,N);for i=2:N-1xx=(i-1)*(X/(N-1);u(1,i)=sin(pi

28、*xx);end;for k=1:Mu(k,1)=0;u(k,N)=0;end;A=zeros(N-1,N-1);%构个M行N列的矩阵用于存放时间t和变量%t=0时刻的赋值，边界条件% x=0,x=1处的边界条件%隐格式的矩阵形式中的A矩阵赋值N=21A(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-2xA(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);

29、L=eye(n);y=zeros(n,1);x=zeros(n,1);for i=1:N-1d(i,1)=sin(pi*(i-1)*(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)

30、;for i=2:ny(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:1x(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=x end for k=1:M u(k,1)=0; end; disp(u); x,t=meshgrid(1:N,1:M); surf(x,t,u); xlabel('x'),ylabel('t'

31、),zlabel('u');title('古典隐格式');%此程序得到图是3-7%追赶法结束%更新右边矩阵%每次外循环更换右边矩阵%显示差分法求得的值%将区域划分成网格对每个点赋值再画图图3-7古典隐格式程序结果图(r=2)古聘隐格或图3-8精确数值解、古典隐格式程序结果的Y-Z平面图(r=2)0.01图3-9古典隐格式在取不同网格比时的误差传播结果图图3-10不同时间取值时精确解、与古典隐格式的值对比图（网格比 r=2）红线表示精确解、 蓝色线表示差分格式的解图3-7古典隐格式在r=2的图形与精确解是一致的。图3-8精确数值解、古典隐格式的Y-Z平面图结果可以

32、看出古典隐格式在t=0.2处的值的边界值和精确解还是有误差的。图3-5是r分别取0.245、0.5、0.72、1.125时的误差传播图像,边界位置网格数为5处的误差为0.01得到的，可以看出r取不同的值时都是稳定的；即古典隐格式对任意的网格比稳定性。从图3-10可以看出隐格式随着时间古旺隐格式-误差侍搐（=0.245古闭隐格式.误差传播司.5古此隙格式误差停皤kCi.72古旺隐格It-误羞侍播E 1250 005-二 0.D06-t=0 OSS时括浦解、古茴混捂式的解招纬.对比t=® OS时.帼醐料、古曲隆格式的解曲跳对此0.711，t=0福确解、古曲格式的酥咛线此比n 折III1=

33、0 175时,指确解、古鼓院捂式的解拍姓.对比 0.25 d _0.01 n DQ05-0-EO的增加，差分格式计算的结果和精确结果越来越大；隐格式虽然对任意网格比都是稳定的，但是计算的精确度是它的不足。3.2.3 Richardson 显格式 close all clc T=0.2 X=1.0000 M=41 N=11 ra=(T/(M-1)/(X/(N-1)A2);fprintf('稳定仕系数S=ra 为:n');disp(ra);u=zeros(M,N); for i=2:N-1xx=(i-1)*(X/(N-1);u(1,i)=sin(pi*xx);end; for k=

34、1:M u(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)=-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)

35、;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);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-1

36、u(2,i)=x(i,1) end %通过隐格式已求得第二行的值u(2,i)，下面是Richardson格式的过程for k=2:M-1%时间轴第3行开始到第 M行for i=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);%区域划分进行赋值画图surf(x,t,u); xlabel('x'),ylabel('t'),zlabel('u'); titl

37、e('Richardson格式');%此程序得到的图形是图3-11Richardscn 格式31 X 10o a图3-12精确数值解、Richardson显格式程序结果的 Y-Z平面图(r=0.5)图3-11 Richardson显格式程序结果图(r=0.5)精匍数佰解图3-13 Richardson显格式在取不同网格比时都不稳定的结果图图3-11是Richardson显格式在r=0.5时的程序结果图，可以看出不稳定。图3-12是精确数值解、Richardson显格式程序结果的Y-Z平面图。从图3-13可 以看出Richardson显格式在取不同网格比时都出现不稳定现象，和理论

38、结果相一 致。所以说Richardson显格式绝对不稳定，这种差分格式不能使用。后面有改进 的格式D-F格式。3.2.4 Crank-Nicholson 格式 close all clc T=0.2 X=1.0000 M=41 N=21 ra=(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-1xx=(i-1)*(X/(N-1);u(1,i)=sin(pi*xx);% 边界条件end; for k=1:M

39、 u(k,1)=0; u(k,N)=0;%边界条件end;%Crank-Nicholson格式需要两个矩阵，下面的A、B,参照Crank-Nicholson格式的矩阵形式%下面对A进行填充赋值A=zeros(N-1,N-1);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-1);B(1,1)=1-ra;B(N-1,N-1)=1-ra;B(1,2)=ra/2;B(N-1,N-2)=ra/2

40、;for m=2:N-2B(m,m-1)=ra/2;B(m,m)=1-ra;B(m,m+1)=ra/2;end;%以下是一一追赶法求%下面对B矩阵进行填充赋值%U的上次对角线即为A的上次对角线%按层外围大循环，即时间步长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(i-1,i);U(i,i)=A(i,i)-L(i,i-1)*U(i-1

41、,i);endfor i=1:N-1d(i,1)=sin(pi*(i-1)*(1.0/(N-1);endfor k=1:M-1m=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-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);for i=2:n

42、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:Mu(k,1)=0;end;disp(u); % u的值全部求出，以下画图x,t=meshgrid(1:N,1:M);surf(x,t,u);xlabel('x'),ylabel('t'),zlabel(

43、'u');title(' Crank-Nicholson "式');%此程序得到的图像是图3-14图3-14 Crank-Nicholson格式程序结果图(r=2)0.90.B 70.6=0.5D.4图3-15精确数值解、Crank-Nicholson格式程序结果的 Y-Z平面图（r=2）0.70.50.40.3g0.1t=fl.O5 8 我布祥、ON格式的就曲践浦比t=0O站时，精醺睬勺C擀式的醉曲技对比E1蹈时椅漏解、CJM梏式的解曲坡对比t=0 175时精砺解、C的梧式的解曲城对比0=3$IRQi jJS图3-17不同时间取值时精确解、与C-N格

44、式的解对比图（网格比r=2）红线表示精确解、蓝色线表示差分格式的解图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时的解对比图，计算还是具有误差。3.2.5 Du Fort Frankle 格式close all clc T=0.2 X=1.000

45、0 M=41 N=21 ra=(T/(M-1)/(X/(N-1)A2);fprintf('稳定仕系数 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%其实初始矩阵已经将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)=

46、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,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

47、); for i=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); 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-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=meshg

48、rid(1:N,1:M);surf(x,t,u);xlabel('x'),ylabel('t'),zlabel('u');title('Du Fort Frankle格式');%此程序为 Richardson格式的改进，得到图 3-18图3-18 Du Fort Frankle格式程序结果图(r=2)0 90 00.7UbU 60 40.30.20 1U图3-19精确数值解、Du Fort Frankle格式程序结果的 Y-Z平面图（r=2）图3-20 Du Fort Frankle格式在取不同网格比时的误差传播结果图吗时，抽琳峰

49、.网信式的岷曲缄利比1=DO»I辑齿斟、Mfert的料启摄对比图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 时误差传播图像，不同的网格比误差都不会扩撒。说明这种格式是稳定的，和理论上的结果相一致。图3-21是不同时间取值时精确解、与 D-F格式的解对比， 随时间的增加计算的值和差分得到的值有误差。此种格式虽然是绝对稳定的，但是计算的精度还是有

50、待提升。四、结论及感想从程序得出的结果与精确解的对比来看和理论是上的结果基本一致。比如古典显格式网格比r小丁等丁 0.5（偏微分方程的系数a取1）才稳定，从MATLAB 编程运行的结果也可以看出r小丁等丁 0.5是稳定的，r大丁 0.5不稳定。乂如 Richardson格式理论上是绝对不稳定的，从编程的结果在取不同的网格比 Richardon格式都是不稳定的，理论和结果一致。理论上对不稳定格式进行改进 使其稳定，比如得到的 D-F格式，取不同的格式网格比都是稳定的，很好的验 证了改进的格式的稳定性。本次报告首先推导了一维抛物线型偏微分方程的一般差分格式。 分别是古典 显格式、古典隐格式、Ric

51、hardson显格式、C-N格式进版的Richardson格式D-F 格式。推导中得到各种格式的截断误差， 从理论上分析了各种格式的稳定性。 对 丁不稳定的格式进行修改得到稳定的格式，即 Richardson格式的修改。通过MATLAB编程实现了各种格式的程序实现。用实验的方法来验证理论 结果，能更好的理解各种差分格式的稳定性及操作过程。报告中的程序都是基本程序，误差图与二维x-u的图像（见附录）都是在基本程序的基础上对参数的修 改得到的图像。通过这次报告的完成学到了很多的东西。首先，对编程有了进一步的了解； 尤其是使用MATLAB编程。在这个过程中也遇到了很多的问题，通过查阅资料 并利用网络

52、资源寻求解决办法。其次，对丁差分格式在程序上的实现并不是简单 的书写，需要步步衔接好；比如好几种格式都用到差分格式的矩阵形式，尤其是隐式格式不能直接求解，需要应用追赶法进行求解；编写过程中通过直接求解得 出的结果都不正确，通过追赶法才能得到正确的结果。最后，我们专业是电磁场与微波技术且偏计算，经常遇到偏微分方程的求解；所以这次试验毫无疑问的对 我们理解有限差分法和求解方程提供了很大的帮助。最后，感谢张老师在课堂上的知识的讲解；同时也感谢在完成报告期间对我 提供帮助的同学！误差传播三维图(以显格式(图 3-5)close allclcT=0.2X=1.0M=41N=8u=zeros(M,N);%

53、构造一个 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.01for k=1:Mu(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);%此处为古典显格式endenddisp(u); %显示差分法求得的值x,t=meshgrid(1:N,1:M);%将区域划分成网格对每个点赋值再画图subplot(2,2,1)surf(x,t,u);xlabel('x'),ylabel('t'