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

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

3、6;)flx-x。)f2x0)(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,ki>u(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、，

4、、1-u2.、3、u(x,tk)=u(xi,tk)()i,k(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

5、)2一一h不得)rho(h)(2-8)得到对位置的二阶偏导数-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=fi'hhk1k,Ji.kkkkUi-Ui2Ui1-'2uiUi1=fih最后得到古典显格式的差分格式为k1kkkkU(1-2ra)5r：

6、工Ui1ut1其中r=J，古典显格式的差分格式的截断误差是o«+h2)h2.1.2古典显格式稳定性分析古典显格式(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

7、-4rasin2j-h-120三ra-121结论：当0WraW，时，所以古典显格式是稳定的。22.2.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-1

8、9）写成kkJUj一山OlTuk书-2心+u；Ik2=fih2jkkJ-.kkkkUj-Uj-斤ujL2UjUp=fj最后得到古典隐格式的差分格式为(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

9、,显然112ra-2racosj二hj=1,2,N-1112ra-2racos(j二h)11+2ra(1cosjh)142j二h14rasin2九H<1包成立。结论：对于Vr0,即任意网格比下，古典隐格式是绝对稳定的2.3Richardson格式2.3.1 Richardson格式的推导将1小书和t=ty,代入式(2-3)得u(xi,tk1)=u(xi,tk)(q)i,k(tk1-tk)o(tk1-tk)2)1t(2-25)u(Xi,tkj)=u(xHtk)(-)i,k(tkd-tk)o(tk-tk)2)二t即f-u2u(xi,tk1)=u(xi,tk)(一)i,k-o():t(2-26

10、)；:u2u(xi,tk)=u(x,tk)-(-)i,k-o(.)Ct由此得到可得(2-27)_22f(X,tk)o(h):uu(x,tk1)-u(xi,tk1)2、(>,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=

11、,，古典显格式的差分格式的截断误差是o(d+h2)h2.3.2 Richardson稳定性分析将Richardson显格式(2-31)写成如下矩阵形式u：1=2厂C-2Iu：u：,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+16ra

12、sin>1(包成立)(2-29)(2-30)(2-31)(2-32)(2-33)(2-34)(2-35)(2-36)(2-37)结论：上式对任意的网比都包成立，即Richardson格式是绝对不稳定的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，k2

13、o(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=tf处的一阶偏导数可以用下式表示: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

14、,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)-'2U(Xj,tk1)U(Xji,tk1)U(Xj1,tk)2U(Xj,tk)U(Xj,tk)|=一2'22ILh2h2所以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)

15、+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)(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+h

16、2)。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二)racosj(h'j=(2-54)(1r:)-racosj(h1-2asi2世=212asfnh2使P(H)<1即1-2asi2-1<12rasin2j2h-1(2-55)22j二h2rh2一1-2asin1r2is4n<1ra2-22jih2(2-56)2j二h2j二h1-2rasin12ras

17、in-222j二h2j二h-1-2rasin1-2rasin22(2-57)2j二h2h-2rasin_2rasin222j二h2j二h-rasin-1-rasin22(2-58)结论：Crank-Nicholson格式对任意网格比也是绝对稳定的。5.DuFortFrankle格式(Richardson格式的改进)将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格式

18、是绝对稳定的。因为此格式是Richardson格式的改进格式，因此截断误差还是o(i2+h2)。3.6总结(1)古典显格式古典显格式的差分格式为u尸=(12ra)u：+r«(u)+u：fj截断误差：o(i+h2)o1稳定性：当0£raM1时，古典显格式稳定。2(2)古典隐格式古典隐格式的差分格式为(12ra)u：-:工1u：1u：=u：，.f；截断误差：o(i+h2)o稳定性：任意网格比古典隐格式绝对稳定。(3) Richardson显格式Richardson显格式的差分格式为u/=2ru(u：书2u：+u：)+qk，+2丁fik截断误差：o(T2+h2)。稳定性：任意网格

19、比Richardson格式绝对不稳定。(4) Crank-Nicholson格式Crank-Nicholson格式的差分格式为1 +。)u：*-卷心；+ukJ1)=(1-r«Ju：+学年书+u：,)+tfjk*+f；)j截断误差：o(2h2)。稳定性：Crank-Nicholson格式对任意网格比绝对稳定。DuFortFrankle格式(1+2ra)uik+=2口(u+uik1)+(1-2ra)uk+2wfik截断误差：o("h2)。稳定性：DuFortFrankle格式对任意网格比绝对稳定。三、MATLAB实现五种差分格式对偏微分方程的求解及误差分析3.1精确数值解一.2

20、二u二u.:t一2-Xt,0u(x,0)=sin(二x)u(0,t)=u(1,t)=0上述偏微分方程的精确解是u(x,t)=eTtsin(二x)(0<x<1,t_0)区域取值范围：0ExEl,0<t<0.20用MATLAB对精确解进行编程画出三维图像精确解程序如下:closeallclcx,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');

21、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('稳定性系数S=ra为:n');disp(ra);%显示网格比fori=2:N-1xx=(i-1)*(X/(N-1);u(1,i

22、)=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);%将区域划分成网格对每个点赋值再画图surf(x,t,u);xlabel('x'),ylabel('t'),zlabel(&

23、#39;u');title('古典显格式');%此程序得到的是图3-3图3-3古典显格式程序结果图（r=0.5）图3-5古典显格式在取不同网格比时的误差传播结果图时，情确颦、古曲显格式的解曲线对比t=0125,横礴解.古典显格式的解曲线对比0.3£i03D.2D40.G口田1图3-6不同时间取值时精确解、与古典显格式的值对比图（网格比r=0.5）红线表示精确解、蓝色线表示差分格式的解图3-1、图3-3对比可以看出，精确解和古典显格式（网格比r=0.5）的图形是一致的。图3-4精确数值解、古典显格式的Y-Z平面图结果可以看出古典显格式的边界值和精确解一样。图3-

24、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=0

25、,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,1

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

27、,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&

28、#39;),zlabel('u');title('古典隐格式');%追赶法结束%更新右边矩阵%每次外循环更换右边矩阵古典隐榨式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

29、得到的，可以看出r取不同的值时都是稳定的；即古典隐格式对任意的网格比稳定性。从图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.

30、3Richardson显格式closeallclcT=0.2X=1.0000M=41N=11ra=(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隐格式的程序一样，只是求一行，此处不

31、再做注释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;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

32、);fori=2:ny(i,1)=u(1,i)-L(i,i-1)*y(i-1,1);endx(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%通过隐格式已求得第二行的值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格式e

33、ndenddisp(u);%显示求解的值x,t=meshgrid(1:N,1:M);%区域划分进行赋值画图surf(x,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是

34、精确数值解、Richardson显格式程序结果的Y-Z平面图。从图3-13可以看出Richardson显格式在取不同网格比时都出现不稳定现象，和理论结果相一致。所以说Richardson显格式绝对不稳定，这种差分格式不能使用。后面有改进的格式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

35、-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格式需要两个矩阵，下面的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);%下面

36、对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)=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的上次对角

37、线即为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%按层外围大循环，即时间步长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=

38、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);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

39、('x'),ylabel('t'),zlabel('u');title('Crank-Nicholson林式')；此程序得到的图像是图3-14图3-14Crank-Nicholson格式程序名果图(r=2)th电第时.，眄薛、c制格式的解曲城对比时扑斛.CH格式的第曲找对比图3-17不同时间取值时精确解、与C-N格式的解对比图（网格比r=2）红线表示精确解、蓝色线表示差分格式的解1t聘时,精福解、C用格式的斜曲城对比图3-14是程序运行得到的Crank-Nicholson格式在网格比取r=2的结果，和精确解图像一致。在t=0.2时

40、从图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 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

41、: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)=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

42、);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%追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,

43、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(1:N,1:M);surf(x,t,u);xlabel('x'),ylabel('t'),zlabel('u');title('DuFortFrankle格式')；此程序为Richardson格式的改进，得到图3-1

44、8图3-18DuFortFrankle格式程序Z果图(r=2)图3-19精确数彳1解、DuFortFrankle格式程序结果的Y-Z平面图（r=2）DuForiF幽kle格式误差悔强（=0245DuForiFankle格式误差拷Jfi(MJ5DuFmFenHe格式误差传JfifO.72DuForiFMk岫格式栽整传播图3-20DuFortFrankle格式在取不同网格比时的误差传播结果图图3-21不同时间取值时精确解、与D-F格式的解对比图（网格比r=2）红线表示精确解、蓝色线表示差分格式的解D-F格式是Richardson格式（绝对不稳定）的改进，从图3-18可以看出当r时D-F格式是稳定的

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

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

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

48、疑问的对我们理解有限差分法和求解方程提供了很大的帮助。最后，感谢张老师在课堂上的知识的讲解；同时也感谢在完成报告期间对我提供帮助的同学！误差传播三维图(以显格式(图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');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处的边界条件fo

49、rk=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'),zlabel('u');title('古典显格式-误差传播r=0.245,);clcT=0.2X=1.0M=41N=1

50、1u=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('古典显格式-误差彳播r=0.5');clcT=0.2X=1.0M=41N=13u=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;for